function model_output=gfop_soe_ces_cap_int_go_solver(V_k,D,X,lambda,beta,int_share,va_share,lab_share,gamma,theta,rho,t_k,s_k,K)

%This function solves for the gains from optimal trade and/or industrial policy,
%starting from the initial equilibrium which is assumed to feature no trade
%or industrial policy.

%CES version with elasticity of substitution rho across sectors
%Inclues intermediate goods and capital
%Capital is assumed to be in fixed supply
%scale economies in gross inputs


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Start with some basic expressions

%Function to calculate domestic lambda_hat change in the cost function
    function lambda_hat=lambda_hat_calc(c_hat)
        lambda_hat=c_hat.^(-theta)./(c_hat.^(-theta).*lambda+1-lambda);
    end

%Function to calculate intermediate price index given p_hat
    function intprice_hat=intprice_hat_calc(p_hat)
        p_hat_rep = repmat(p_hat',K,1);
        intprice_hat=prod(p_hat_rep.^int_share,2);
    end

%Function to calculate c_hat given w_hat,l_hat,r_hat,intprice_hat
    function c_hat=c_hat_calc(w_hat,l_hat,r_hat,intprice_hat)
        c_hat_num=(((w_hat.^lab_share).*r_hat.^(1-lab_share)).^va_share).*intprice_hat;
        c_hat_denom=(1+s_k).*((l_hat.*(w_hat).^(1-lab_share.*va_share).*(r_hat).^((lab_share-1).*va_share).*intprice_hat.^(-1)).^gamma);
        c_hat = c_hat_num./c_hat_denom;
    end

%Function to calculate b_hat (change in beta) given p_hat
    function b_hat=b_hat_calc(p_hat)
        b_hat=p_hat.^(1-rho)/(sum(p_hat.^(1-rho).*beta));
    end

%Function to calculate p_hat given c_hat
    function p_hat=p_hat_calc(c_hat)
        p_hat=(c_hat.^(-theta).*lambda+1-lambda).^(-1./theta);
    end

%Calculations leading to total revenue

%Function to calculate subsidy revenue (first term)
    function subsidy_revenue=subsidy_revenue_calc(lambda_hat,E_prime)
        subsidy_revenue=-sum(s_k.*lambda.*lambda_hat.*E_prime); %Note this is negative
    end
%Function to calculate export tax revenue (net of subsidy costs)
    function tariff_revenue=tariff_revenue_calc(c_hat)
        tariff_revenue=sum((t_k.*(1+s_k)-s_k).*(1-t_k).^(theta).*c_hat.^(-theta).*X);
    end
        
%Function to calculate R_prime
    function R_prime = R_prime_calc(subsidy_revenue,tariff_revenue)
        R_prime=subsidy_revenue+tariff_revenue;
    end
        
%Calculations leading to total expenditure

%Function to calculate RHS vector
    function rhs_exp=rhs_exp_calc(w_hat,r_hat,c_hat,b_hat,tariff_revenue)
        a_exp=b_hat.*beta*(w_hat*sum(lab_share.*V_k)+r_hat*(1-sum(V_k.*lab_share))+D);
        b_exp = b_hat.*beta*tariff_revenue+int_share'*((1+s_k).*(1-t_k).^(theta+1).*c_hat.^(-theta).*X);
        rhs_exp=a_exp+b_exp;
    end

%Function to calculate matrix for inversion
    function lhs_matrix_exp=lhs_matrix_exp_calc(lambda_hat,b_hat)
        a1=-s_k.*lambda_hat.*lambda;
        A1=repmat(a1',K,1);
        b1=b_hat.*beta;
        B1=repmat(b1,1,K);
        AB=A1.*B1;
        C=(int_share.*repmat((1+s_k).*lambda_hat.*lambda,1,K))';
        lhs_matrix_exp=eye(K)-(AB+C);
    end

%Function to calculate change in expenditure
    function E_prime=E_prime_calc(rhs_exp,lhs_matrix_exp)
        E_prime=lhs_matrix_exp\rhs_exp;
    end
       
%Compute initial total expenditure from the data
E=(V_k./va_share-X)./lambda;

%Compute expenditure an alternative way
    function E_prime_other=E_prime_other_calc(w_hat,r_hat,R_prime,E_prime,c_hat,b_hat)
        aa=b_hat.*beta*(w_hat*sum(lab_share.*V_k)+r_hat*(1-sum(V_k.*lab_share))+D+R_prime);
        bb=int_share'*(lambda.*lambda_hat.*E_prime+(1-t_k).^(theta+1).*c_hat.^(-theta).*X);
        E_prime_other=aa+bb;
    end

%Calculation of l_hat given the rest of the functions

    function l_hat=lhat_calc(w_hat,lambda_hat,E_prime,c_hat)
        %lhat1=(va_share.*(1+s_k))./(V_k*w_hat);
        %lhat2=lambda_hat.*lambda.*E_prime;
        %lhat3=(1-t_k).^(theta+1).*c_hat.^(-theta).*X;
        l_hat=(va_share.*(1+s_k))./(V_k*w_hat).*(lambda_hat.*lambda.*E_prime+(1-t_k).^(theta+1).*c_hat.^(-theta).*X);
    end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solution Algorithm

%Initial values
w_hat=1;
r_hat=1;
l_hat=ones(K,1);
p_hat=ones(K,1);
intprice_hat=ones(K,1);
b_hat=ones(K,1);

%Pick initial c_hat to reflect subsidies
c_hat=c_hat_calc(w_hat,l_hat,r_hat,intprice_hat);

%Max Iterations
max_c_hat_iter=15; % 
max_w_iter=5000;    %also iterations for r_hat

%Error tolerance
max_e_w=.00001;
max_e_c_hat=.0001;

%Step sizes
ss_w=.05;
ss_c_hat=.05;


%%%%%%%%%%%%%
%Algorithm
for i=1:max_w_iter
    for j=1:max_c_hat_iter
            %Start by computing implied endogenous variables
            p_hat=p_hat_calc(c_hat);
            lambda_hat=lambda_hat_calc(c_hat);
            b_hat=b_hat_calc(p_hat);
            tariff_revenue=tariff_revenue_calc(c_hat);
            rhs_exp=rhs_exp_calc(w_hat,r_hat,c_hat,b_hat,tariff_revenue);
            lhs_matrix_exp=lhs_matrix_exp_calc(lambda_hat,b_hat);
            E_prime=E_prime_calc(rhs_exp,lhs_matrix_exp);
            subsidy_revenue=subsidy_revenue_calc(lambda_hat,E_prime);
            R_prime = R_prime_calc(subsidy_revenue,tariff_revenue);
            l_hat=lhat_calc(w_hat,lambda_hat,E_prime,c_hat);
            E_prime_other=E_prime_other_calc(w_hat,r_hat,R_prime,E_prime,c_hat,b_hat);
            %Avoiding any negative values of l_hat
            %l_hat= max(l_hat,.001*ones(K,1));
            %Now compute new implied c_hat
            intprice_hat=intprice_hat_calc(p_hat);
            c_hat_new=c_hat_calc(w_hat,l_hat,r_hat,intprice_hat);
            %Check the error
            error_c_hat=c_hat_new-c_hat;
            Z_c_hat=sum(abs(error_c_hat));
              if Z_c_hat<max_e_c_hat
                break;
              end
            %Update c_hat by a step
            c_hat=c_hat+ss_c_hat*(c_hat_new-c_hat);
    end     
       %Update wages based on excess demand for labor
       excess_lab_dem=sum(l_hat.*V_k.*lab_share)-sum(V_k.*lab_share);
       Z_w = abs(excess_lab_dem);
       w_hat=w_hat+ss_w*(w_hat*excess_lab_dem);
       iterw=i;
       %Update interest rate
       r_hat = w_hat/(1-sum(V_k.*lab_share))*sum(l_hat.*V_k.*(1-lab_share));
       if Z_w<max_e_w
           break;
       end                 
end  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Model Output

%Calculating the change in the price index
P_hat=sum(p_hat.^(1-rho).*beta)^(1/(1-rho));

%Hat change in welfare
welfare=(w_hat*sum(V_k.*lab_share)+r_hat*(1-sum(V_k.*lab_share))+D+R_prime)/((1+D)*P_hat);


%New trade shares
new_lambda=lambda.*lambda_hat;



model_output.new_lambda=new_lambda
model_output.lambda_hat=lambda_hat;
model_output.p_hat=p_hat;
model_output.c_hat=c_hat;
model_output.l_hat=l_hat;
model_output.intprice_hat=intprice_hat;
model_output.b_hat=b_hat;
model_output.tariff_revenue=tariff_revenue;
model_output.rhs_exp=rhs_exp;
model_output.lhs_matrix_exp=lhs_matrix_exp;
model_output.E_prime=E_prime;
model_output.E=E;
model_output.R_prime=R_prime;
model_output.subsidy_revenue=subsidy_revenue;
model_output.w_hat=w_hat;
model_output.r_hat=r_hat;
model_output.Z_w=Z_w;
model_output.E_prime_other=E_prime_other;
model_output.welfare=welfare;
model_output.P_hat=P_hat;

%Model output for Harberger Analysis
va_new = w_hat*sum(l_hat.*V_k);    %total value added at the new equilibrium
go_new = w_hat*l_hat.*V_k./va_share;   %gross output by sector at the new equilibrium
model_output.GO_VA_share=go_new./va_new;    %ratio of gross output by sector to total value added at new equilibrium

z_hat = (l_hat.*(w_hat).^(1-lab_share.*va_share).*(r_hat).^((lab_share-1).*va_share).*intprice_hat.^(-1));   %hat change in gross input
model_output.GO_real = 1-1./z_hat;  %gross input reallocation

model_output.fs_share=(1-t_k).^(theta).*c_hat.^(-theta).*X/va_new;   %foreign expenditures on domestic goods, as a share of domestic value added

fz_hat = (1-t_k).^(theta+1).*c_hat.^(-theta -1);    %hat change in effective gross inputs supplied to foreigners
model_output.fGO_reallocation=1-1./fz_hat;  % (negative) percent change in effective gross inputs supplied to foreigners

end